# Load all snow depth and density data
# Load snow pit data
load("all_pits")
# Load Camera Data
load("all_cam")
# AWS depth
load("all_hr_AWS")
# JW SNOTEL Date
load("JW_snotel")
# Import SPICE Data
spice <- read.csv("AWS/SPICE/BACKUP.CSV") %>%
select(c(2,3)) %>%
mutate(date = sub("T.*", "", Timestamp), time = str_remove(sub(".*T", "", Timestamp), "Z"),
Timestamp= format(as.POSIXct(paste(date, time), tz = "MST"), "%Y-%m-%d %H:%M:%S"),
date = date(Timestamp) ,site = "UF_uc",
depth = Height..m. - 2.47) %>%
filter(depth <= 2.5 & depth >= 0) %>%
mutate(depth = replace(depth, depth <= 0.01, 0)) %>%
group_by(date) %>%
summarize(sd_sonic = median(depth, na.rm = TRUE)*100) %>%
filter(date >= "2021-10-15") %>%
mutate(rolling_sd = rollmean(sd_sonic, 3, fill = NA), site = "UF_UC")
#Plot camera snow depths
SD_Daily <- ggplot(all_cam) +
geom_line(aes(x = date, y = depth_avg, color = site ))+
labs(title = "Camera Snow Depth", x = "Date", y = "Snow Depth (cm)")
ggplotly(SD_Daily) %>% layout(hovermode = "x unified")
#Plot spice snow depths
spice_sd <- ggplot(spice) +
geom_line(aes(x = date, y = sd_sonic, color = "SD"))+
geom_line(aes(x = date, y = rolling_sd, color = "3 day rolling SD"))+
labs(title = "Spice Snow Depth", x = "Date", y = "Snow Depth (cm)")
ggplotly(spice_sd) %>% layout(hovermode = "x unified")
# Mutate data to include SD, Density, and/or SWE
# Snow Pit Density
pit_density <- all_pits %>%
mutate(density_avg = rowMeans(.[,5:7], na.rm = TRUE), weight = densityTop - densityBot) %>%
group_by(date, site) %>%
summarise(bulk_density = weighted.mean(density_avg, weight, na.rm = TRUE), pit_SD = max(heightAG, na.rm = TRUE))
## `summarise()` has grouped output by 'date'. You can override using the
## `.groups` argument.
# Daily Wx sonic snow depth data, merge with correct density values
daily_sonic <- all_hr %>%
select(c(1,2,10)) %>%
mutate(date = date(TIMESTAMP), site = ifelse(site == "bf", "BF_Wx", "UF_Wx")) %>%
group_by(date, site) %>%
summarise(depth = mean(DBTCDT_Avg, na.rm = TRUE)*100) %>%
merge(., pit_density, all = TRUE) %>%
filter(site == "UF_Wx" | site == "BF_Wx") %>%
select(-c(pit_SD)) %>%
mutate(site = ifelse(site == "BF_Wx", "BF_Wx_sonic", "UF_Wx_sonic"))
## `summarise()` has grouped output by 'date'. You can override using the
## `.groups` argument.
# Merge daily SPICE data with correct density values
spice_density <- merge(spice, pit_density, all = TRUE) %>%
mutate(site = replace(site, site == "UF_UC", "UF_Wx")) %>%
filter(site == "UF_Wx") %>%
group_by(date, site) %>%
summarise_all(sum, na.rm = TRUE) %>%
select(-c(pit_SD)) %>%
mutate(site = "UF_UC", bulk_density = replace(bulk_density, bulk_density == 0, NA))
# Daily Camera data
daily_cam <- all_cam %>%
select(-c(1:4)) %>%
relocate("date", .before = "site")
# Daily Joe Wright SWE Data
JW_SWE <- JW_snotel %>%
select(c(12, 19, 13)) %>%
rename(SWE = snow_water_equivalent)
#Plot Snow pit density
density <- ggplot(pit_density) +
geom_line(aes(x = date, y = bulk_density, color = site ))+
labs(title = "Snow Pit Density", x = "Date", y = "Density (kg/m^3)")
ggplotly(density) %>% layout(hovermode = "x unified")
# Merge all automated data with snow pit depth and density
AWS_daily <- list(daily_cam, pit_density, daily_sonic, JW_SWE, spice_density) %>%
reduce(full_join, by = c("date", "site")) %>%
mutate(sd = coalesce(depth_avg, pit_SD, depth, rolling_sd),
sd = ifelse((site == "UF_Wx_sonic" & date >= "2022-06-14") |
(site == "BF_Wx_sonic" & date >= "2022-06-03") |
(site == "BF_Wx" & date == "2022-06-03"), 0, sd),
bulk_density = coalesce(bulk_density.x, bulk_density.y, bulk_density)) %>% #coalesce to fill gaps in camera and sonic time series with pit
group_by(site) %>%
mutate(density = na.approx(bulk_density, rule =2, na.rm = FALSE), .after = sd, SWE = ifelse(is.na(SWE), density*sd/100, SWE),
burn = ifelse(startsWith(site, "BF_"), "yes", "no")) %>%
# extrapolate of density calculate SWE
select("date", "site", "SWE", "sd", "density", "burn") %>%
relocate("burn", .after = "site") %>%
filter(date < "2022-06-15")
# Plot the timeseries data
#SD
SD_Daily <- ggplot(AWS_daily) +
geom_line(aes(x = date, y = sd, color = site )) +
labs(title = "Automated SD", x = "Date", y = "Snow Depth (cm)")
ggplotly(SD_Daily) %>% layout(hovermode = "x unified")
#SWE
SWE_Daily <- ggplot(AWS_daily) +
geom_line(aes(x = date, y = SWE, color = site )) +
labs(title = "Automated SWE", x = "Date", y = "SWE (mm)")
ggplotly(SWE_Daily) %>% layout(hovermode = "x unified")
#SWE
Density_Daily <- ggplot(AWS_daily) +
geom_line(aes(x = date, y = density, color = site )) +
labs(title = "Density Timelapse", x = "Date", y = "SWE (mm)")
ggplotly(Density_Daily) %>% layout(hovermode = "x unified")
AWS_daily_summary <- AWS_daily %>%
group_by(site) %>%
summarize(max = max(SWE, na.rm = T))
# SWE increases
snow_change <- AWS_daily %>%
arrange(date) %>%
group_by(site) %>%
mutate(accum = ifelse(SWE > lag(SWE), SWE - lag(SWE), NA),
ablat = ifelse(SWE < lag(SWE), lag(SWE) - SWE, NA))
# accumulation
cam_accum <- ggplot(snow_change, aes(x = date, y = accum, fill = site)) +
geom_bar(stat = "identity") +
facet_wrap(~site, ncol = 2) +
labs(title = "New Snow Accumulation", x = "Date", y = "Accumulation (mm)") +
theme(legend.position = "none")
ggplotly(cam_accum) %>% layout(hovermode = "x unified")
## Warning: Removed 1107 rows containing missing values (position_stack).
#ggsave("newsnowaccum.png")
# ablation
cam_ablat <- ggplot(snow_change, aes(x = date, y = ablat, fill = site)) +
geom_bar(stat = "identity") +
facet_wrap(~site, ncol = 2) +
labs(title = "Snow Ablation", x = "Date", y = "Ablation (mm)") +
theme(legend.position = "none")
ggplotly(cam_ablat) %>% layout(hovermode = "x unified")
## Warning: Removed 946 rows containing missing values (position_stack).
#ggsave("snowablation.png")